home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / MacMETH 3.2.3 / More Examples / Newton.MOD < prev    next >
Text File  |  1995-12-13  |  6KB  |  240 lines
  1. MODULE Newton; (* HS 19.10.91 *)
  2.  
  3.   FROM SYSTEM IMPORT Exp, Ln, Sqrt, Sin, Cos, ArcTan;
  4.   FROM Terminal IMPORT BusyRead;
  5.   FROM CursorMouse IMPORT GetMouse;
  6.   FROM InOut IMPORT Read, Write, WriteString, WriteLn;
  7.   FROM Windows IMPORT SetWindow, ResetWindow;
  8.   FROM GraphicWindows IMPORT Window, OpenGraphicWindow, CloseGraphicWindow, Clear,
  9.                              SetPen, MoveTo, Dot, IdentifyPos;
  10.  
  11.  
  12.   CONST NrPixels = 256;
  13.  
  14.         (* colours *)
  15.         white   = 449;
  16.         yellow  = 65;
  17.         green   = 321;
  18.         zyan    = 257;
  19.         blue    = 385;
  20.         magenta = 129;
  21.         red     = 193;
  22.         black   = 1;
  23.  
  24.  
  25.   VAR ch             : CHAR;
  26.       v,w            : Window;
  27.       environ        : BITSET;
  28.       ox,oy,nx,ny    : INTEGER;
  29.       dx,dy,max,p,q  : REAL;
  30.       xmax,xmin      : REAL;
  31.       ymax,ymin      : REAL;
  32.       ux,uy,uux,uuy  : REAL;
  33.       deltax         : INTEGER;
  34.       maxiter        : INTEGER;
  35.       this           : INTEGER;
  36.  
  37.  
  38.   PROCEDURE f(x: REAL): REAL;
  39.   BEGIN
  40.     RETURN x * Cos(6.0*x) * Exp(-x*x)
  41.   END f;
  42.  
  43.   PROCEDURE ForeColour(c : LONGINT); CODE 0A862H;
  44.   PROCEDURE BackColour(c : LONGINT); CODE 0A863H;
  45.  
  46.   PROCEDURE Colour(k : INTEGER);
  47.   VAR c : LONGINT;
  48.   BEGIN
  49.     IF k >= maxiter THEN c := black
  50.     ELSE
  51.       k := k MOD 16;
  52.       CASE k OF
  53.           0 : c := white;
  54.         | 1 : c := yellow;
  55.         | 2 : c := yellow;
  56.         | 3 : c := green;
  57.         | 4 : c := green;
  58.         | 5 : c := zyan;
  59.         | 6 : c := zyan;
  60.         | 7 : c := blue;
  61.         | 8 : c := blue;
  62.         | 9 : c := blue;
  63.         |10 : c := magenta;
  64.         |11 : c := magenta;
  65.         |12 : c := magenta;
  66.         |13 : c := red;
  67.         |14 : c := red;
  68.         |15 : c := red;
  69.       END;
  70.     END;
  71.     ForeColour(c);
  72.   END Colour;
  73.  
  74.   PROCEDURE GetPos(VAR i,j: INTEGER);
  75.   CONST ML = 15;
  76.   VAR mouse : BITSET; x,y : INTEGER;
  77.   BEGIN mouse := {};
  78.     REPEAT GetMouse(mouse,x,y) UNTIL NOT(ML IN mouse);
  79.     REPEAT GetMouse(mouse,x,y) UNTIL ML IN mouse;
  80.     REPEAT GetMouse(mouse,x,y) UNTIL NOT(ML IN mouse);
  81.     IdentifyPos(w,x,y);
  82.     i := x; j := y;
  83.   END GetPos;
  84.  
  85.   PROCEDURE SANE(VAR e: BITSET; OpWord: CARDINAL); CODE 0A9EBH;
  86.  
  87.   PROCEDURE SaveFPEnv;
  88.   BEGIN SANE(environ, 3)
  89.   END SaveFPEnv;
  90.  
  91.   PROCEDURE ClearFPEnv;
  92.   VAR e: BITSET;
  93.   BEGIN e := {}; SANE(e, 1)
  94.   END ClearFPEnv;
  95.  
  96.   PROCEDURE RestoreFPEnv;
  97.   VAR e: BITSET;
  98.   BEGIN e := environ; SANE(e, 1)
  99.   END RestoreFPEnv;
  100.  
  101.  
  102.   PROCEDURE Calculate;
  103.   VAR x,y: REAL; ix,iy: INTEGER;
  104.  
  105.     PROCEDURE iteration(pe,qe,xe,ye: REAL) : INTEGER;
  106.     VAR a,y1,y2,x1,x2: REAL; counter: INTEGER;
  107.  
  108.     BEGIN
  109.       IF xe = ye THEN RETURN maxiter END;
  110.       x1 := xe;
  111.       x2 := ye;
  112.       counter := 0;
  113.       y1 := f(x1);
  114.       y2 := f(x2);
  115.       REPEAT
  116.         a := y1 - y2;
  117.         IF a = 0.0 THEN RETURN maxiter END;
  118.         a := a / (x1 - x2);
  119.         x2 := x1; y2 := y1;
  120.         x1 := x1 - y1/a;
  121.         y1 := f(x1);
  122.         INC(counter);
  123.       UNTIL (ABS(y1) < max) OR (counter >= maxiter);
  124.       RETURN counter
  125.     END iteration;
  126.  
  127.   BEGIN (* Calculate *)
  128.     SaveFPEnv;
  129.     ClearFPEnv;
  130.     dx := (xmax-xmin) / FLOAT(NrPixels);
  131.     dy := (ymax-ymin) / FLOAT(NrPixels);
  132.     y := ymin;
  133.     FOR iy := 0 TO NrPixels-1 DO
  134.       x := xmin;
  135.       FOR ix := 0 TO NrPixels-1 DO
  136.         this := iteration(p,q,x,y);
  137.         SetWindow(w);
  138.         Colour(this);
  139.         Dot(w,ix,iy);
  140.         ResetWindow;
  141.         x := x + dx;
  142.       END;
  143.       y := y + dy;
  144.       BusyRead(ch); IF (ch = ' ') THEN iy := NrPixels END;
  145.     END;
  146.     RestoreFPEnv;
  147.   END Calculate;
  148.  
  149.  
  150.   PROCEDURE ShowFunction(v: Window);
  151.   VAR ix,iy: INTEGER; x,y,y0,fmin,fmax,ry: REAL;
  152.       a: ARRAY [0..NrPixels] OF REAL;
  153.   BEGIN
  154.     SaveFPEnv;
  155.     ClearFPEnv;
  156.     fmin := 0.0; fmax := 0.0;
  157.     FOR ix := 0 TO NrPixels DO
  158.       x := xmin + FLOAT(ix) * ( (xmax-xmin) / FLOAT(NrPixels) );
  159.       y := f(x); a[ix] := y;
  160.       IF y > fmax THEN fmax := y
  161.       ELSIF y < fmin THEN fmin := y
  162.       END;
  163.     END;
  164.     FOR ix := 0 TO NrPixels DO a[ix] := a[ix] - fmin END;
  165.     FOR iy := 0 TO NrPixels DO
  166.       SetWindow(v);
  167.       ForeColour(zyan);
  168.       Dot(v, NrPixels DIV 2, iy);
  169.       ResetWindow;
  170.     END;
  171.     ry := ABS(fmax-fmin);
  172.     y0 := ABS(fmin) / ry * FLOAT(NrPixels);
  173.     iy := TRUNC(y0);
  174.     FOR ix := 0 TO NrPixels DO
  175.       SetWindow(v);
  176.       ForeColour(zyan);
  177.       Dot(v, ix, iy);
  178.       ResetWindow;
  179.     END;
  180.     FOR ix := 0 TO NrPixels DO
  181.       y := (a[ix] / ry) * FLOAT(NrPixels);
  182.       iy := TRUNC(y);
  183.       SetWindow(v);
  184.       ForeColour(red);
  185.       Dot(v, ix, iy);
  186.       ResetWindow;
  187.     END;
  188.     RestoreFPEnv;
  189.   END ShowFunction;
  190.  
  191.  
  192. BEGIN
  193.  
  194.   max := 1.0E-5;
  195.   maxiter := 255;
  196.   xmax := 4.0; xmin := -4.0;
  197.   ymax := 4.0; ymin := -4.0;
  198.   p := 0.0; q := 0.0;
  199.  
  200.   OpenGraphicWindow(v,40,20,NrPixels+4,NrPixels+20,"Funktion",Clear);
  201.   Clear(v);
  202.   ShowFunction(v);
  203.  
  204.   OpenGraphicWindow(w,330,100,NrPixels+2,NrPixels+20,"Newton",Clear);
  205.   Clear(w);
  206.  
  207.   LOOP
  208.     WriteString ('drawing Newton.'); WriteLn;
  209.     Calculate; ch := 0C;
  210.     WriteString ('zoom requested Y/N:');
  211.     Read(ch); Write(ch); WriteLn;
  212.     IF CAP(ch) # "Y" THEN EXIT END;
  213.     WriteString ('define window !'); WriteLn;
  214.     GetPos(ox,oy);
  215.     ux := xmin + FLOAT(ox)*dx;
  216.     uy := ymin + FLOAT(oy)*dy;
  217.     GetPos(nx,ny);
  218.     deltax := ABS(nx-ox);
  219.     uux := xmin + FLOAT(nx)*dx;
  220.     uuy := uy + FLOAT(deltax)*dy;
  221.     xmin := ux; xmax := uux;
  222.     ymin := uy; ymax := uuy;
  223.  
  224.     SetWindow(w);
  225.     ForeColour(black);
  226.     SetPen(w,ox,oy);        MoveTo(w,nx,oy);
  227.     SetPen(w,nx,oy);        MoveTo(w,nx,oy+deltax);
  228.     SetPen(w,nx,oy+deltax); MoveTo(w,ox,oy+deltax);
  229.     SetPen(w,ox,oy+deltax); MoveTo(w,ox,oy);
  230.     ResetWindow;
  231.  
  232.     GetPos(ox,oy); (* wait for mouse click *)
  233.     Clear(w);
  234.   END;
  235.  
  236.   CloseGraphicWindow(w);
  237.   CloseGraphicWindow(v);
  238.  
  239. END Newton.
  240.